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Abstract 

The extended framework of Hamilton's principle and the mixed convolved action 
principle provide new rigorous weak variational formalism for a broad range of initial 
boundary value problems in mathematical physics and mechanics. In this paper, their 
potential when adopting temporally higher order approximations is investigated. The 
classical single-degree-of-freedom dynamical systems are primarily considered to 
validate and to investigate the performance of the numerical algorithms developed 
from both formulations. For the undamped system, all the algorithms are symplectic 
and unconditionally stable with respect to the time step. For the damped system, 
they are shown to be accurate with good convergence characteristics. 
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Introduction 

Despite of its origin in particle dynamics, Hamilton's principle (Hamilton 1834; Hamilton 
1835) has been with us for a long time throughout broad range of mathematical physics 
(Bretherton 1970; Gossick 1967; Landau & Lifshitz 1975; Slawinski 2003; Tiersten 1967). 
However, it suffers from two main difficulties such as (i) use of end-point constraints and 
(ii) adoption of Rayleigh's dissipation for non-conservative systems. The first difficulty 
relates to the proper use of initial conditions resulting from the restrictions on the func- 
tion variations. In Hamilton's principle, the variations vanish at the end points of the time 
interval, which, in turn, implies that the functions are known at these two instants. For a 
typical dynamic problem, one does not know how the considered system evolves at the 
end of the time interval. Usually, this is the main objective of the analysis, which means 
that there may be a serious philosophical or mathematical inconsistency in Hamilton's 
principle. Second difficulty relates to the inability to incorporate irreversible phenomena. 
Hamilton's principle itself only applies to conservative systems. With Rayleigh's 
dissipation (Rayleigh 1877), irreversible processes can be brought into the framework of 
Hamilton's principle. However, this approach is not satisfactory in a strict mathematical 
sense, since the variation of Rayleigh's dissipation enters in an ad-hoc manner. 

Historically, to resolve such difficulties in Hamilton's principle, Tonti (Tonti 1973) 
suggested that convolution should replace the inner product for variational methods in 
initial value problems. Somewhat earlier, (Gurtin 1963; Gurtin 1964a,b) introduced the 
convolution functional, and could reduce the initial value problem to an equivalent 
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boundary value problem. However, the functional by Gurtin is complicated and it never 
can recover the original strong form. Following the ideas of Tonti and Gurtin, 
Oden and Reddy (Oden and Reddy 1976) extended the formulation to a large class 
of initial boundary problems in mechanics, especially for Hellinger-Reissner type 
mixed principles. 

More recently, Riewer (Riewe 1996; Riewe 1997) adopted the use of fractional 
calculus to accommodate dissipative dynamical systems. This is an attractive idea, and 
many other researches including Agrawal (Agrawal 2001; Agrawal 2002; Agrawal 2008), 
(Atanackovic et al. 2008), (Baleanu and Muslih 2005), Dreisigmmeyer and Young 2003; 
Dreisigmeyer and Young 2004), (El-Nabulsi and Torres 2008), and (Abreu and Godinho 
2011) have proposed similar approaches. However, surprisingly, none of these papers 
include an analytical description validating their approach for the most fundamental 
case, a classical Kelvin-Voigt single-degree-of-freedom (SDOF) damped oscillator. 

Recently, two new variational frameworks for elastodynamics such as extended 
framework of Hamilton's principle (EHP, (Kim et al. 2013)) and mixed convolved action 
principle (MCAP, (Dargush and Kim 2012)) were established by using mixed variables. 
While EHP adopts a mixed Lagrangian formalism given in (Apostolakis and Dargush 
2012; Apostolakis and Dargush 2013a; Sivaselvan et al. 2009; Sivaselvan and Reinhorn 
2006), it provides a new and simple framework that correctly accounts for initial 
conditions within Hamilton's principle. EHP resides in an incomplete variational frame- 
work since it requires Rayleigh's function for dissipative systems and cannot define the 
functional action, explicitly. On the other hand, MCAP clearly resolves long-standing 
problems in Hamilton's principle. With MCAP, a single scalar functional action 
provides the governing differential equations, along with all the pertinent boundary and 
initial conditions for conservative and non-conservative linear systems. Thus, in 
theoretical aspects, MCAP is certainly preferred rather than EHP, however, there still 
remains a challenge for MCAP to have the generalized framework of other than linear 
problems. While EHP can be numerically implemented for viscoplasticity continuum 
dynamics, MCAP currently suffers from finding the explicit functional action for that 
problem. Since both methods provide sound basis to develop various space-time finite 
element methods for linear initial boundary value problems, here, the focus is initially 
on investigating their potential when employing higher-order temporal approximations. 

The remainder of the paper is organized as follows. Next, in Section New variational 
formalisms, some relevant background on EHP and MCAP are provided, especially for 
the SDOF Kelvin-Voigt system. In Section Numerical implementation, discretization 
scheme and numerical algorithms are provided when temporally higher-order approxi- 
mations are adopted in both approaches. Basic numerical properties of the developed 
methods are closely examined in Section Basic numerical properties. Then, some 
numerical examples are presented to investigate and to validate all of these developed 
algorithms for practical problems of the forced vibration in Section Numerical examples. 
Finally, some conclusions are provided in Section Conclusions. 

New variational formalisms 

In this section, new variational frameworks for the SDOF Kelvin-Voigt system displayed 
in Figure 1 were reviewed for the development of higher order temporal finite element 
methods from both approaches. 
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Figure 1 SDOF Kelvin-Viogt damped oscillator. 



With mass m, damping coefficient c, the known applied force f[t) with time t, and 
stiffness k = 1/a with a representing the flexibility, EHP and MCAP could formulate 
the variational framework for this model in terms of the displacement of the mass u(t) 
and the impulse of the internal force /(f) in the spring. 

Weak form for the Kelvin-Voigt model in EHP 

Following the ideas in (Sivaselvan and Reinhorn 2006), the EHP associated with this 
problem defines Lagrangian L and Rayleigh's dissipation <p as 

L(u,u,j;t)=-mii 2 + -aj 2 -ju+fu (1) 



and 

<p(u;t) = 1 -c[u(t)} 2 (2) 

where a superposed dot represents a derivative with respect to time. 
Then, the functional action A for the fixed time interval from t 0 to t is given by 

t 

A(u,u,j;t) = J L(u,u,J;t) dr (3) 

to 



and, in EHP, the first variation of A is newly defined as 

t t 

8 Anew = -8 j L(u,u, J;r)dr + j ^ .' — -r5w<ir+ \mu r5«j =0 (4) 

to t 0 

by adding the counterparts (the underlined terms in Eq. (4)) to the terms without the 
end-point constraints in Hamilton's principle. 

Such added terms have effect on confining a dynamical system to evolve uniquely 
from start to end with the unspecified values at the ends of the time interval such as ii 
(to), w(to), u(t) and u(t). Then, interpreting the unspecified initial terms as sequen- 
tially assigning the known initial values completes this formulation. Thus, in EHP, the 



Kim SpringerPlus 2014, 3:458 
http://www.springerplus.eom/content/3/1/458 



Page 4 of 25 



given initial velocity Uq is assigned first, and the given initial displacement Uo is 
assigned next by 

u(t Q ) = Uq (5) 

and 

Su(t a ) = 0 ( ■:u(t 0 ) = u 0 ) (6) 

The subsequent zero-valued term (6) needs not appear explicitly in the new action 
variation, so that the new definition (4) with the sequential assigning process such as 
(5) and (6) can properly account for the initial value problems. It should be noted that 
in EHP, the dependent initial condition J 0 can be identified by 

mu Q + cu 0 +/ 0 -/ 0 = 0 (7) 

where j 0 is the initial internal impulse of the known applied force/ given by 

to 

Jo = J f(r) dr (8) 

In Eq. (8), the time interval [-°°, t 0 ] is used to represent that this is the time interval 
before the initial time we are considering. 
To check this, let us substitute Eqs. (1), (2) into Eq. (4). Then, we have 

£ 

8 Anew = - J \m u 8u-c ii 8u-j 8u +f 8u + a J 8J -u 5/j dr + \mu 5wj = 0 

to 

(9) 

Doing integration by parts on muSu, aj 8J , and -u8J in Eq. (9) yields 

8A NEW = \mu(t) 8u(t) -mii(to) 8u(to) \-[m u Su] 1 ^ + [(u-aj) 8j\ t 



j [m u + cu + j -f^jSu dr + j (aj- u)8J dr = 0 



(10) 



For arbitrary variations of 8u and <5/for the time interval [t 0 , t], the governing differential 
equations are given by 

m u + cu + J -f = 0; a j-u = 0 (11) 

along with constitutive relation as 

u-aj = 0 (12) 

With the underlined terms in Eq. (10), the trajectory of the damped oscillator is 
firstly uniquely confined by 

k{t) =u{t); 8u{t) =8u(t) (13) 

while the given initial conditions are identified sequentially by Eq. (5), Eq. (6) and Eq. (7). 

Thus, with EHP, Hamilton's principle can account for compatible initial conditions to 
the strong form. It is not a complete variational method, since it still requires the 
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Rayleigh's dissipation for a non-conservative process and the first variation of the 
functional action cannot yield the proper weak form explicitly. However, the framework 
is quite simple and it can be readily applied to problems other than linear elasticity 
with the use of Rayleigh's dissipation. 

For a representative example, let us consider SDOF elasto-viscoplastic model in 
Figure 2. Rayleigh's dissipation to define rate-deformation for the slider-dashpot k\ can 
be given by 



?>(/; t) 



(14) 



in terms of Macaulay bracket ( • ) and absolute value of / whereby y\ and F y represent 
viscosity and yield force, respectively. Thus, in EHP, the action variation for this model 
is defined by adding up Eq. (4) and SA lp 



SA„, 



a/ 



-SJdT 



J\(\j\- F y) sga{j)SJdT 



(15) 



and 



(16) 



where the underlined term represents the rate-deformation for the slider-dashpot, 

Note that the adding terms (16) are the counterparts to the terms without the 
end-point constraints in Hamilton's principle that obtained from the compatibility 
condition 



a] + U\ = u 



(17) 



With Eq. (4) and Eqs. (14, 15 and 16), the governing differential equations for 
Figure 2. 



mu + cii+j=f; aj-ii + ui 



0 



(18) 



are properly recovered in EHP along with proper initial conditions such as Eqs. (5, 6 and 7) 
and «i at t 0 . 




Kim SpringerPlus 2014, 3:458 
http://www.springerplus.eom/content/3/1/458 



Page 6 of 25 



Weak form for the Kelvin-Voigt model in MCAP 

As well described in (Dargush and Kim 2012), MCAP defines the convolved action for 
the SDOF Kelvin-Voigt damped oscillator as 

A^u, u, it, J, J ,/; tj =~m(u*u)-~a(j */) + *u^j + —c(u*u)-(u* f^-u(t)j(0) 

(19) 

where a superimposed arc represents a temporal left Riemann-Liouville semi-derivative. 
Referred to (Oldham and Spanier 1974; Samko et al. 1993), this is defined by 

t 

u{t) = fi3^„V t ) 3 _i_l / J^dt (20) 
w V 0 ) KI 1(1/2) dt J {t-rf 2 1 ; 



where T( ) denotes the Gamma function. 

In Eq. (19), the symbol * represents the convolution of two functions over time, such 
that 

t 

(<p*<f>)(t)= I <p(T)4>(t-T)dr (21) 



Meanwhile, the last term /(0) in Eq. (19) represents the initial impulse corresponding 
to f(t) that is given by 



/(*) = / fir) dr (22) 

In MCAP, the stationarity of the action (19) yields the following weak form in time 

SA = m{Su *u)-a(Sf * /) + US J *u^j + (^Su* J J + c(Su *u)-(Su * fJ-Su(t)j(0) = 0 

(23) 

After performing classical and fractional integration by parts on the appropriate 
terms in Eq. (23) as follows (Apostolakis and Dargush 2012), we have 

SA = (du* {mii + cii +/-/}) + {SJ * {-aj + ii}) + Su{t) |raii(0) + c»(0) + 7(0) -](0)} 
+ Su(0){mu(t)} + dJ(t){-a/(0)+u(0)}-SJ(0){-aj(t)} = 0 

(24) 

For the sake of completeness, the fractional integration by parts formula is given 



t t 



[Df v ) (t) (D){U) (f - t) dr = |g (t) <P(t - r) dr + g>(0) 4>{t) (25) 



o o 
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For arbitrary variations of u and /, Eq. (24) emanates the governing differential 
equations in mixed forms as 

m u + cu + J -f = 0; a j-ii = 0 (26) 

along with the proper initial conditions 

wm(0) + cw(0)+/(0)-/(0) = 0; -«/(0) + w(0) = 0 (27) 

Note that the initial variations such as Su{Q) and (5/(0) vanish due to Eq. (27). In other 
words, in MCAP, we can identify the dependent initial conditions such as /(0) and /(0) 
from the usual given initial conditions w(0) and w(0) as well as the known initial 
impulse ;'(0). 

As shown in Eq. (24) and Eqs. (26, 27), every governing equations and initial 
conditions are satisfied weakly in MCAP, where it incorporates both conservative 
and non-conservative components within the unified functional action (19). Thus, 
it resolves the long-standing problem in Hamilton's principle. However, MCAP still 
requires a generalized framework to embrace various irreversible phenomena. In 
particular, currently, it does not have the functional action for the problem shown in 
Figure 2. Also, it should be noted that any pair of complementary order of fractional 
derivatives in Eq. (19) yields Eqs. (26, 27) due to the integration by parts property of 
complementary order of fractional derivatives 

t t 
J (t) (£& 4>){t-r)dT= J^(t) 4>{t-r) dr + <p(0) 0(f) (28) 

0 0 



for 0 < a < 1. 



Numerical implementation 

The weak form (9) in EHP and the weak form (23) in MCAP include, at most, first 
derivatives of the primary variables u{t) and J{t) as well as the variations Su{t) and SJ{t). 
Consequently, we have C° temporal continuity requirement on primary variables and 
the variations, thus, there are many cases to develop higher order temporal finite 
element methods. As we shall see in this section, three kinds of quadratic temporal 
finite element methods in each framework are developed, since they are practically 
sufficient and accurate in computational aspects as discussed next. The numerical 
methods developed here are summarized in Table 1. 

Table 1 Developed quadratic temporal finite element methods in each framework 

Algorithms Description 

Jquad J{t) and 6J(t): quadratically approximated. 

u(f) and 8u{t) : linearly approximated. 
Uquad u(f) and 6u{t): quadratically approximated. 

/(f) and 5J(t): linearly approximated. 
UJquad u(f) and 6u(t)\ quadratically approximated. 

J{f) and 5J(t): quadratically approximated. 
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Algorithms from EHP 

By introducing the fixed time step h for each time duration, that is, t r = r h, Eq. (9) can 
be written 



N 



SA = Yj SA >- 

= V^- I / \m u 8ii-c ii Su-j Su+f Su + a J SJ-u 8j\dT-\p Suf t ' ' i I =0 



V,-i 



(29) 



where SA r represents the action variation in the r* time duration [t r _i, t r ]. Also, 
p represents linear momentum, where 

p(t) = m u{t) (30) 

For t r _ i < r < i„ temporally linear shape functions such as L r _ 1 at t r _ ! and Z.,. at i r 
are given by 

I r _ 1 (r)=i(i r -r) (31) 

L r (r) =i(r-^_i) (32) 



Also, by introducing the center point t c for the time interval [t r _i, f r ] as 
(f r -f r _i) // 



(33) 



2 2 

temporally quadratic shape functions Q r _ i at t r _ i, Q r at t„ and Q c at t c can be written as 

2 
h 

2 



Qr-iW = 7j(r-*r)(r-t c ) (34) 
Q r {r)=^{j-t r ^){r-t c ) (35) 



Q e M=-72(T-fr)(T-*r-l) ( 36 ) 

h 

With linear temporal shape functions (31)-(32) and quadratic temporal shape 
functions (34)-(36), we can develop every algorithms of EHP presented in Table 1. 
For a representative case, Jquad algorithm can be obtained from the main approxima- 



tions as 




7(T) = Q r -i(T)/,._i + Q r (T)Jr + Q c {t)Jc 


(37) 


SJ{t) = Q^tfSJr-i + Q r (T)SJ r + Q c (t)SJ c 


(38) 


u(t) = L r _i(r) u r -\ + L r (r) u r 


(39) 


8u{t) = L r _i(r) Su r -i + L r (r) Su r 


(40) 


f(T)=L r - 1 (r)f r _ 1 +L r (T)f r 


(41) 


and the subsequent approximations as 




/(T) = Qr-l(r)7r-l + Qr(r)Jr + Q c {?)Jc 


(42) 
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SJ(t) = Q r _ 1 (T)dJ r . 1 + Q r (j)SJ r + Q c (r)SJ c 
u(t) = L r _i(r) u r -i + L r {r)u r 
8u{t) = L r - i(r) Su r -i + L r {r) Su r 



(43) 
(44) 
(45) 



Substituting Eqs. (37, 38, 39, 40, 41, 42, 43, 44 and 45) into Eq. (29), and integrating 
yields 



SA r 



yyi c i 5 2 fa " fa * \ 

■ — {u r -U r -l}+-{u r -U r -l} + Pr-^r-\+-^Jr--Jc--^fr-l-^fr\ Su >- 



f 1 7 8 1 1 

/ 8 8 16 \ 2 2 

\ " 3F^ " 3T /r + 3l 7 7 + 3 ^ " 3 Mr 



-1 + g«r )<S/r 



5/ c = 0 



(46) 



By making the coefficient of (<5w r _i, r5w„ SJ r -\, 8]„ SJ C ) equal to zero in Eq. (46), we 
have four independent equations given by 

m c 5 1 2 h " h * \ 

— {K r - M r _i} + - {U r - Ur-x} -p r _ 1 --/ r -i + ~J r + ~J C ~ ~f r _ x - -f r J =0 (47) 



yyi c 1 5 2 fa ^ fa ^ 

— {lir-Ur-i} +-{u r -U r -i} + p r - -J r -1 + ^Jr ~ ^Jc ~ r -i " 3/r ) = 0 ( 48 ) 



3h Jr - 1 + 3h Jr 3h Jc 



1 7 8 

Sh 1 "- 1 + 3h Jr ~ 3h Jc 



6 6 



1 5 

— U r -i -\- — Uj 



0 



(49) 
(50) 



While deriving Eqs. (47, 48, 49 and 50), the equation from the underlined term 
in (46) is discarded because it is not independent, which can be obtained from 
adding Eq. (49) and Eq. (50). 

From either Eq. (49) or Eq. (50), we can express J c in terms of J r _ 1, /„ 1 and u r Then, 
replacing J c in the other independent equations with the equation of JJJ r -\,J„ 
yields the matrix equation of Jquad algorithm as 



1 (X + 6cha) 

12 ha 

1 (-X + 6cha) 

12 ha 



2 
1 

2 

2a 



1 (X + 6cha) . 1 

12 ha 2 
1 {-X+6cha) 1 

12 ha 2 
2a 

0 -T 




(51) 



where X is given by 



X= Ylma-h 2 



(52) 
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Similarly, we have the Uquad algorithm as 

l(3m + ch) 1 (Y + 6cha) 

6 h 2 
1 (-Y + 6cha) 

~6 " 



3 h 
l(-9m + 2ch) 

3 h 

2 (6 m + ch) 

3 h 
l{9m + 2ch) 

3 h 
1 (-3m + ch) 

~3 h 

2 (-6m + ch) 

3 hT 



h 2 
2 X 

IF 

1 (Y + 6cha) 

6 h 2 

1 (-Y+6cha) 

6 

2 X 
"3 /? 





I 3 



6 fr-l 
h 

fr-1 + o/r 



(53) 



where Y is given by 
r = 24ma + h 2 
Also, we have the UJquad algorithm as 



1 (I2ma + 4:cha + h 2 ) q 1 (Y + 6cha) 



12 fea 
1 (-36wa + 8c/2« + /2 2 

12 /za 

2 (6m + ch) 

3 A 

1 (36wa + 8c/2a-/z 2 ) 

12 ha 
1 (-12ma + 4c/za-/z 2 

12 ha 

2 (-6m + ch) 

3 hT 



6 h 2 
1 (Y-6cha) 

6 P 

2 ho 

3 A 2 

i (y + 6c/?a) 

6 " ~h 2 
1 (Y-6cha) 

6 P 

2 HQ 

3 /z 2 




(54) 




(55) 



with the adequate substitution of u c and J c in terms of J r _ 1( /„ w r _ 1( and « r . 



Algorithms from MCAP 

Previously, MCAP was numerically implemented through linear temporal shape 
functions for classical SDOF oscillators and systems that utilize fractional-derivative 
constitutive models by (Dargush 2012). Here, continuing through this line, but, the 
quadratic temporal finite element methods are developed. 

As well described in (Dargush 2012), for any non-negative integer m and n, we have 
the following relation 



for the convolution of the semi-derivatives of power functions. 

To evaluate the convolution of semi-derivatives of polynomial shape functions, here, 
Eq. (56) is frequently used. 
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Since we cannot have summation form of the action variation in convolution integral 

N 

(that is, 8A&y ] SA r ), let us consider the action variation over one time-step [0, h] as 

r=l 

SA(u, u, u, J, J , /; t—>}i j = 

m(Su *u)-a (SJ */) + ( SJ *u^j + (Su */ J + c (5 2" * 2") - ^<5m *fj - Su(t)j(0) = 0 

(57) 

where temporally linear and quadratic shape functions of t (0 < t < h) are defined as 



L 0 {t) = l- t l 


(58) 




(59) 


2 ht+ 2) 


(60) 




(61) 


Q c (t) = -^(t 2 -ht) 


(62) 



Then, subsequent approximations are given by 



Lo(t)=-l 


(63) 




(64) 




(65) 




(66) 


Q c (t) = -^(2t-h) 


(67) 



Now, let us consider Jquad algorithm for a representative one. 

With approximations (58)-(67), the convolution component (sj *u) in Eq. (57) can 
be written as 



(sj *«)(*) = [SJ 0 Sh $Jc\ 



L<57o Sh SJ C \ 



Qi*£o Qi*£i 
Q c *l 0 Q c *Ii 



(68) 



in terms of row vector |_ • J, matrix [•], and column vector {•}. 
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Each component of matrix in Eq. (68) can be directly evaluated by using Eq. (56). For 
a representative one, [Q 0 * L 0 ) (t) is computed as 



(Qo*A))(f) 



-ht- 



* D, 



(t) 



2 



1 f 3 



3t 2 h 1 h 



(69) 



-^-~ht + - — + — -~t 



h 3 2 



2 2 2 2 



Then, by letting f— »/z in Eq. (69) due to the underlined term in Eq. (57), Eq. (69) 
yields 



(Qo*I o )(^0 



2 



9 1 t 3 3, 3 i 2 h 2 h 

t-- ht-\ 1 1 

h 3 2 22 2 2 



2 



1A 3 3., 3 A 2 /z 2 h 



h 3 2 



-/z/z- 



2 2 2 2 



(70) 



Following the same procedures as in Eqs. (69, 70), one finds 



8J*u= [8J Q Sh 8Jc\ 



1 1 

~6 6 

5 1 

6 6 

2 2 

~3 3 



(71) 



In a similar way, 



S u* J = [Suq Su\\ 



15 2 



6 6 
1 1 



6 6 



/o 



(72) 



and for the viscous dissipation term 



c(5 u* u) = [8uo 8ui\ 



2 2 
c c 



MO 
Ml 



(73) 



With evaluation of typical integer order convolution components in Eq. (57), we have 
the following discretized weak form of Jquad: 
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[Suq 8u\ J 



- m 


m- 


h 


h 


YYl 


YYl 


-~~h 


~h - 




- 1 




~6 


8J C \ 


5 




6 




2 




-~3 


- c 


c-i 



\8} 0 8J, 8J C \ 



a 


la 


8a 


~3h 


~3h 


3h 


la 


a 


8a 


~3h 


~3h 


3h 


8a 


8a 


16a 


3h 


3h 





Jo 

Jc 







]l 


>!:) 




UJ 



+ [8uq Sui J 



2 2 
c c 



h h 

6 3 

/z h 

3 6 



1 5 

~6 6 3 
112 



(74) 



With the known initial conditions u 0 and To, the variations Su 0 and <5/ 0 vanish. Thus, 
the weak form reduces to the following: 



m 


m 


f «0 1 
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«o 



8u\ 



112 

6 6 3. 



8a' 
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(75) 



8u\ 



h h 
L3 6 



|°}-*«i[/„] =0 



Then, grouping the terms according to the variations and allowing the arbitrary 
variations on 8u\, 8J\, 8J C , one obtains following equations 



YH C II 

■ — {u 0 - Ml} + - {U 0 + Mi} + -Jo + -/i 

n L 6 6 



^-(7/o+/i-8/ c )+^M 0 + ^i =0 
3« 6 6 

^(-8/o-8/i + 16/ c )-^ Mo + ^«i = 0 



2 J- h f - h f 

3 Jc 3/0 g/i 



; 0 



(76) 
(77) 
(78) 



Again, with the adoption of the same strategy as Eqs. (47, 48, 49, 50 and 51) in EHP 
to express J c in terms of / r _ i, /„ « r _ 1( and w„ finally, we have 



1 (X+6cha) 1 
12 /zfl 2 



1 (X-6c/z«) 

12 /zfl 
1 



o 



(79) 



where X is defined in Eq. (52) and Qj is given by 



n - n - 
Qa = 3/0 + ^/1+70 



(80) 
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More generally, for the n th time step with t n = nh, one may write the Jquad algorithm 
of MCAP 
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12 ha 
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a 


~2 
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U n -i 
Jn 



, , + l o 



(81) 



where 



/z - /z - , 
Qj„ = 3-^"-i + +/«-! 



(82) 



Similarly, we can develop the Uquad algorithm as 



1 (3m + c/z) 1 {Y + 6cha) 

~3 h 6 

2 (6m + c/z) 2 (X) 

3 A ~3 ~W 

1 (9m + 2c/z) 1 (r-6/z 2 + 6c/z 

3 h 6 " 

2 (-6m + c/z) 

3 jT~ 



u„ 1 
/„ J 



/z 2 
2 (X) 

'3 /z 2 



/m-1 



* ; ft. 



(83) 



where X and I' are given respectively in Eq. (52) and Eq. (54), while Q„ is given by 

Qu„=P„-i+L-i (84) 
Also, we have the UJquad algorithm as 



1 (mX + 6c/zam + c 2 /z 2 a) 1 (6m + ch) 

12 ham 12 m 

1 (6m + c/z) 1 (X) 

12 m 12 hm 

1 (mX-6cham+ c 2 h 2 a) 1 (-6m + c/z) 

12 liam 12 



ham 
1 (-6m + c/z) 

12 m 



m 
12 /zm 



+ 



/z 2 - h 2 f 



(85) 



where 

fh c/z 2 \~ /ft cd ! \ ; , 
^=(3 + 24mK 1 + (6 + 24mj^ + ^ ^ 



Basic numerical properties 

For the SDOF Kelvin-Voigt model, every algorithm from EHP and MCAP can be 
written in matrix form as 
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Table 2 Algorithms from EHP for the conservative system 
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A 1 x n =A 0 x n _i+f n 

or simply 

x n =A D x n _i+Aj 1 f n 



(87) 



(88) 



where 



A D =A 1 1 A 0 



(89) 



Symplectic nature 

For the undamped case with no external forcing (conservative harmonic oscillator), 
Eqs. (87, 88 and 89) reduce to 



Aleft *n — A r ight x n - 1 

x n =Ax n _i 

A = A-left A rig ht 

where A left and A right in each algorithm are identified in Table 2 and Table 3. 
Table 3 Algorithms from MCAP for the conservative system 
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Table 4 Eigenvalues of A in EHP algorithms 



Algorithms Eigenvalues 



A, = 1 

Jquad 6 m a-2 h 2 ± / \/36 h 2 m 0-3 h 4 

A" = ; ri 

6ma + h 

A, = 1 
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A, = 1 _ 
UJquad h 4 -60amh 2 + 144m 2 a 2 ± i(12h)\]2am-h 2 \^am 

13 ~ h 4 + \2amh 2 + 144m 2 o 2 



In each Table, X and 1^ are given respectively in Eq. (52) and Eq. (54), while Z is 
given by 

h 2 



6ma 



(93) 



Notice that every algorithm shown in Table 2 and Table 3 is time reversible. One can 
exactiy recover the state n - 1 from the state n by setting h—>-h,n—>n-l, and n - 1 — > n. 
For the representative one, one can obtain Uquad algorithm in MCAP as 
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(96) 



which is the exactly same as the Uquad algorithm given in Table 3. 
While deriving Eq. (96), the following relation is used 

-Y + 4X = Y-6h 2 



(97) 



The stability and dissipative character of each developed method can be determined 
by considering the eigenvalues of A in Eq. (92), and the eigenvalues of each method are 
presented in Table 4 and Table 5, respectively. 

Table 5 Eigenvalues of A in MCAP algorithms 
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Eigenvalues 
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Figure 3 Period elongation property of each method. 



Notice that the magnitude of all the eigenvalues including complex conjugate pairs in 
each Table is exactly equal to 1, which can be written simply as 



1 



(98) 



Consequently, in addition to being time reversible, all the presented quadratic 
temporal finite element algorithms are also symplectic, energy conserving, and un- 
conditionally stable for the undamped case. 



Period elongation property in each method 

To check the period elongation property in each developed method, the method by 
(Bathe 1996; Bathe and Wilson 1972) is used for free vibration of the undamped oscillator, 
where the ratio of the time-step h to the natural period T n is a control parameter. Also, 
Newmark's constant average acceleration method and Newmark's linear acceleration 
method are adopted for the references. 

As shown in Figure 3, the numerical dispersion property from EHP and MCAP is 
exactly the same as Newmark's linear acceleration method, when either the primary 
variable u or / is quadratically approximated. On the other hand, when u and / are 
quadratically approximated, UJquad algorithm in each method has the same numerical 
dispersion property better than Newmark's linear acceleration method. Note that all 



Table 6 Numerical simulation cases 
Sinusoidal loading f(t) = 100 sin (10 f) 



El-Centro loading 



J" (i) h = 0.10 
I (ii) h = 0.05 
[ (iii) h = 0.01 

while damping coefficient c = 0.2/r is fixed to deliver a 
non-dimensional damping ratio £ =0.05. 



(i) f = 0.05 

(ii) ( = 0.03 

(iii) { = 0.01 

while the time step is fixed as h = 0.02. 
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Analytical 
h-0.01 
h-0.05 
h=0.10 




0.5 I 1.5 2 2.5 3 3.5 4 4.5 
Time ( 

Figure 4 Displacement history results from Newmark's linear acceleration method. 



the developed methods are unconditionally stable, while Newmark's linear acceleration 
method is a conditionally stable algorithm with the criterion h/T„ < 0.551 

In computational aspects, compared to Newmark's constant average acceleration 
and Newmark's linear acceleration method, all the developed computational 
methods seem practically sufficient and accurate, since they have symplectic, 
unconditionally stable, and less or equivalent period elongation properties, and this 
is the main reason that only quadratic temporal finite element methods are 
developed here. 

Numerical examples 

For all of the numerical examples considered here, with no loss of generality, the 
model parameters are taken in non-dimensional form. In particular, let m = 1 and 



Analytical 
h=0.01 
h=0.05 
h-0.10 




0.5 I 1.5 2 2.5 
Time i 

Figure 5 Displacement history results from Jquad algorithm in EHP. 
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-4< 1 ' ' 1 ' 1 ' 1 

0 0.5 I I.S 2 2.5 3 3.5 4 4,5 

Time i 

Figure 6 Displacement history results from Uquad algorithm in EHP. 



a = l/(4n 2 ), thus, providing a natural period T n = 1 in the SDOF Kelvin-Voigt 
damped oscillator. 

Two loading cases with zero initial conditions are considered for numerical simula- 
tion. The first one is an applied force in the form/(i) = f 0 sin(o> 0 i) with f 0 = 100 and 
<u 0 = 10, and the other is 1940 El-Centro loading. The additional parameters for each 
loading case are summarized in Table 6. 

For the references, the results obtained from each developed method are 
compared to an exact solution for the sinusoidal loading, while the results from 
Newmark's linear acceleration method in OpenSees (Mckeena et al. 2013; McKenna 
2011) are additionally provided. For El-Centro loading, the results from each 
developed method are compared to those from Newmark's linear acceleration 
method in OpenSees. 
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Analytical 
h-0.01 
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h (I.I!) 
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Time t 

Figure 8 Displacement history results from Jquad algorithm in MCAP. 



Simulation results under sinusoidal loading 

Figure 4 displays the numerical solution of displacement versus time, based upon New- 
mark's linear acceleration method, while Figures 5, 6, 7, 8, 9 and 10 are obtained from 
the developed algorithms. 

As seen from the results, all the developed methods have better convergence charac- 
teristics compared to Newmark's linear acceleration methods under sinusoidal loading. 
In particular, UJquad algorithm in each framework shows the most accurate results. 

Simulation results under 1940 El-Centro loading 

The results from 1940 El-Centro loading analysis are displayed in Figures 11, 12, 13, 
14, 15 and 16. In each figure, the Uquad and Jquad algorithms yield the exactly same 
results, while there are slight differences between the newly developed methods and 



Analytical 
h-0.01 
h-0.05 
h=O.I0 




Figure 9 Displacement history results from Uquad algorithm in MCAP. 
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Figure 10 Displacement history results from UJquad algorithm in MCAP. 



Newmark's linear acceleration method. In practical aspects, these differences seem 
negligible, but, note that all the developed methods are unconditionally stable that it 
may be advantageous to have the outlined results before the detailed analysis with the 
new methods. 

Conclusions 

In recent papers, through mixed formulation, two new variational frameworks such as 
EHP and MCAP were formulated for dynamical systems. Theoretically, MCAP is 
preferred to EHP, because unlike previous variational approaches, MCAP does not 
require any dissipation function with ad-hoc rules for taking variations, restrictions on 
the variations at the ends of the time interval, and external specification of initial 
conditions. However, there still remains a challenge for MCAP to have a generalized 



Newmark linear ace 

J quad 

Uquad 

UJquad 
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Time I 

Figure 11 Results from EHP algorithms for El-Centro loading analysis (1% damping ratio) 
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framework embracing various irreversible phenomena. On the other hand, EHP has a 
relatively simple framework: the action variation is newly defined by adding the 
counterparts to the terms without the end-point constraints in Hamilton's principle, 
which confines a dynamical system to evolve uniquely from start to end. Interpreting 
these additional terms as sequentially assigning the known initial values completes this 
formulation. It should be noted that EHP is not a complete variational method, since it 
still requires the Rayleigh's dissipation for a non-conservative process and it cannot 
define the functional action explicitly. Since both mixed formalism provide a rigorous 
foundation to develop various temporal finite element methods for linear elasticity, in 
this paper, their potential when adopting temporally higher order approximations is 
investigated for the classical SDOF Kelvin-Voigt damped system. 




-0.015 1 1 1 1 1 1 1 1 ' 

0 5 III 15 20 25 30 35 40 

Time I 

Figure 13 Results from EHP algorithms for El-Centro loading analysis (3% damping ratio). 
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Figure 14 Results from MCAP algorithms for El-Centro loading analysis (3% damping ratio). 



With the consideration of computational aspects, three quadratic temporal finite 
element methods are essentially developed from each mixed formalism. All the 
developed methods are symplectic and unconditionally stable for the undamped 
conservative harmonic oscillator. Also, from period elongation property studies, it is 
checked that all the developed methods are equivalent or superior to Newmark's linear 
acceleration method that is conditionally stable. For damped forced vibrations, all the 
developed methods are shown to be robust and to be accurate with good convergence 
characteristics. It should be noted that since the new methods utilize mixed formula- 
tions, there exists an inherent disadvantage in a significant increase of the degrees of 
freedom against Newmark's methods when dealing with other than SDOF systems. 
However, this may be somewhat compensated by the general characteristics of a mixed 
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Figure 16 Results from MCAP algorithms for El-Centro loading analysis (5% damping ratio). 



formulation and its broad applicability (Casciaro and Cascini 1982; Commend et al. 
2004; Glowinski et al. 1989; Lee and Filippou 2009). 

As the original Hamilton's principle has been adopted in various applications, the 
applicability of EHP and MCAP are quite broad, spanning many fields of mathematical 
physics and engineering. Future work will be directed toward development of a general- 
ized framework of MCAP, and applications of both formalisms to various engineering 
problems, following the ideas in (Fried 1969; Hulbert 1992; Hulbert and Hughes 1990; 
Li and Wiberg 1996; Pitarresi and Manolis 1991; Bar-Yoseph 1989; Apostolakis and 
Dargush 2013b). 
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